Collision-induced friction in the motion of a single particle on a bumpy inclined line 
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By means of Molecular Dynamics simulations, we investigate the elementary process of avalanches 
and size segregation by surface flow in 2 dimensions: a single ball confined to moving along an 
inclined line consisting of balls. The global characteristics of the motion depend strongly on the size 
of the moving ball relative to the size of the balls on the line, as well as the distribution of the balls 
on the line. We find that in the steady state the friction force acting on the ball is independent of 
material properties like the Coulomb friction coefficient and the coefficient of restitution. Contrary 
to previous notions about the details of the motion, we find that it is very regular and consists of 
many small bounces on each ball on the line. As a result of this regularity, introducing a random 
spacing between the balls on the line has mainly the same influence as a regular spacing of adequate 
length. The insensitivity of the steady state velocity to material properties and to the detailed 
arrangement of the balls on the line allows for an analytical estimation of the mean velocity which 
fits the simulation results very well. We find that results from the 2D case can probably not be 
transferred to the 3D case of a ball moving on a rough inclined plane as easily as has been suggested 
previously. 
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I. INTRODUCTION 

The flow of granular materials has been studied exten- 
sively both experimentally and theoretically, due to its 
ubiquity in nature and its industrial importance. Never- 
theless, many properties of granular flow are still poorly 
understood. Here, we deal with a special case of gran- 
ular flow, namely that along an inclined rough surface. 
Examples are flow in inclined chutes M and all kinds of 
avalanche processes, such as rockslides |g[ , which involve 
segregation processes (known as "inverse grading" to ge- 
ologists), an important phenomenon often encountered 
in granular materials [||. In rotating drums, avalanche 
processes are the motor of size segregation [Q-Q . Segre- 
gation in surface flow is strongly related to the question 
of stability of granular flow on rough surfaces, i.e. to the 
determination of the limiting conditions for the existence 
of a steady state where the flow is neither stopped nor 
accelerated. The rougher the surface encountered by the 
moving particles, the slower they flow and may even come 
to a stop. Thus, the larger the flowing particles are, com- 
pared to the roughness or bumpyness of the surface on 
which they flow, the farther they should travel and accu- 
mulate (in the case of flow on a sand heap) at the bottom 
of the slope. 

The threshold for the onset of granular shear flows as 
well as the limiting conditions for which a stable flow 
(i.e. a steady state) still exists are insufficiently under- 
stood and hard to determine in realistic situations ||. 
In particular the hysterctic properties of granular mate- 
rials, manifesting themselves for example in the differ- 
ence between the static and dynamical angle of repose 



of sandpiles ||, complicate the situation. Though the 
static threshold for the onset of motion on a rough plane 
(i.e. the tilt of the surface large enough to set a resting 
mass of granular material in motion) has been quite thor- 
oughly investigated in experiments and computer simu- 
lations pO| , pl| , the dynamical case, i.e. the stability of 
the flow, is still poorly understood. 

In order to investigate the dynamical situation, we fol- 
low the course taken up by Riguidel et al. pd|- p^| in their 
experimental studies. We consider the elementary pro- 
cess of an avalanche: a single ball of radius R moving 
down an inclined plane onto which other balls of radius r 
are glued. This fixes the roughness of the surface. Such a 
system has been investigated experimentally and numer- 
ically in 3D (i.e. for a ball moving on a plane) |pd,p^,^4| 
and in 2D (for a ball moving on a line of balls) [ 13^ 15|| . 
The simulations in 3D were restricted to the determina- 
tion of the static angle of stability for a ball resting on 
the plane, whereas all experiments and the simulations 
on the line also dealt with the dynamical situation of 
a ball that starts on the plane with some initial veloc- 
ity, which is the situation we are interested in. In all 
cases, three types of motion could be observed, depend- 
ing on the size ratio $ = R/r and on the inclination 
angle of the plane. These results could be summed up 
in a "phase diagram" (typical for both plane and line), for 
which we will later give an example from our simulations. 
The three types of motion observable in experiments and 
simulations are characterized in the following way. In 
regime A, the ball gets trapped on the plane, indepen- 
dently of the initial velocity with which it is launched 
onto it. In regime B, it reaches a constant average ve- 
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locity v in the direction along the plane, in regime C it 
accelerates throughout the whole length (2 m in the ex- 
periments) of the plane, accompanied by visible jumps. 
In the constant velocity regime, the mean velocity v was 
found to be proportional to <I> Q sm9 in 3D, with an expo- 
nent a w 1.3 |oJ1aH , whereas in 2D v 2 depended lin- 
early on sinO ]15| , |16| . In 2D, no simple power law could 
be found for the ^-dependence of the mean velocity. The 
relation v 2 ~ sin 9 was already derived using very gen- 
eral assumptions by Bagnold for flow of many particles 
on an inclined plane |l7j], but no assumption about the 
dimensionality of the system was made there. 

An obvious difference between the 2D and 3D case is 
the fact that in 3D, the particle moving down the plane 
will be deflected in the direction perpendicular to the 
plane inclination, either by rolling down the crooked val- 
leys formed by the balls on the plane, or by obliquely 
impacting a sphere on the plane. We will come back to 
the importance of this possibility, which is absent from 
the simulations we present here, after having discussed 
the details of the motion of the ball moving down the 
line. Recently, a 2D stochastic model has been proposed 
which reproduces the angle dependence of the velocity in 
3D, but gives a different exponent a Jig] ]. 

Here, we discuss the mechanism by which the ball 
keeps a constant velocity on the inclined line and how, 
from the understanding of this mechanism, the transition 
to the stopping and accelerating regime can be explained. 
We will show that the 2D case is generically different from 
the 3D case, i.e. that the effect of disorder on the plane 
cannot be modeled by disorder on a line. Possible reasons 
for the success and the adequateness of the 2D stochastic 
model of ref. |Q in describing the 3D case are discussed. 

The outline of the paper is as follows. After presenting 
our simulation method in section II, we will first show 
in section III that the simulations reproduce experimen- 
tally observed macroscopic behaviour. We then proceed 
to a detailed analysis of the microscopic properties of the 
motion and give an explanation of the mechanism stabi- 
lizing the motion of the ball in the 2-dimensional case. 
The influence of material properties on the motion is in- 
vestigated. On the basis of our simulation results, we 
present in section IV a simplified model for the motion 
of the ball which allows the analytical derivation of the 
mean velocity in the 2-dimensional well as the 

determination of the boundary between regions A and B 
in the phase diagram. The results of this model agree 
with the simulation results. Possible reasons for the dif- 
ference between the 2-dimensional and the 3-dimensional 
system are discussed. 



II. SIMULATION METHOD 

We model the 2D case by the Molecular Dynamics 
(MD) technique p9[ , which was first introduced to the 
simulation of granular materials by Cundall and Strack 



p(| . The MD technique consists of time-integrating New- 
ton's equations of motion for a system of grains starting 
from a given initial configuration. Since our simulations 
are 2-dimensional, the grains have only three degrees of 
freedom, two translational, one rotational. Two grains 
of radii positions Pi, velocities Vi, and angular ve- 
locities u>i (i = 1,2), are in contact when their (virtual) 
overlap £ = max(0, R\ + R2 — \f*2 — P{\) is larger than 
zero ( "soft" grains) . Two unit vectors ft and s are used 
to decompose the forces and velocities into normal and 
shear components: 



r2-n 

n = -p zri 

((n) y , - (n) x J . 
The forces between grains are then given by 

Fij F. n fi -\- F s s 1 

where 



F, 



min(|7 s u s |, \(iF n \) ■ sign(u s ). 



(1) 
(2) 

(3) 

(4) 
(5) 



Here /j, denotes the Coulomb friction coefficient. The rel- 
ative normal velocity v n and the relative shear velocity 
v s (i.e. the relative velocity of the surfaces at the point 
of contact) are defined as 



v n = {v 2 -Vi)-H 

v s = [v% - vi) ■ s + ujiRi + w 2 i?2- 



(G) 
(7) 



A number of different force laws is commonly used in 
MD simulations of granular materials; a discussion of 
their properties can be found in Our choice of F n 

corresponds to a simple linear spring-dashpot, the tan- 
gential force F s is the Coulomb friction law for sliding 
friction, which was regularized for small v s to avoid the 
discontinuity of the Coulomb law at v s = 0. The tan- 
gential damping constant 7 S should have a sufficiently 
high value such that the case F s = — j s v s occurs only for 
very small v s . Only then the interpretation of F s a mere 
regularization of Coulomb friction holds. This force law 
has the advantage of holding equally well in free impacts 
of spheres [2l|j22]| and in long lasting contacts [^3| if the 
particles can be considered as rough hard spheres (i.e. if 
tangential elasticity can be neglected). It also gives a 
velocity dependent coefficient of tangential restitution, a 
feature found to be important in a stochastic model of 
the situation (HHJ] . 

Throughout our simulations we used the parameters 
k n = 2 • 10 6 N/m, 7s = 100 kg/s, [i = 0.13, r = 5 mm 
and M = ^irR 3 p for the mass of the rolling ball with 
p = 7.8 g/cm 3 . The values of these parameters were cho- 
sen to match the steel balls used in |l3|,[l6| ; the choice of 
k n leads to a collision time of the order of 10~ 5 s, which 
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is a typical value for steel balls of this size. The value 
of 7 S is high enough such that F s is reasonably close to 
the exact Coulomb friction law in the sense explained 
above. The damping 7 n is determined by fixing the nor- 
mal coefficient of restitution e„ = — u^/u^, defined by the 
ratio of final and initial normal velocities. Unless stated 
otherwise, the results presented in the next section were 
obtained using e„ = 0.7 and the above mentioned param- 
eters. The influence of the material parameters e n and 
fi on the behaviour of the system will be discussed later. 
The only external force acting on the ball is gravity, the 
gravitational acceleration in the x (y) direction is given 
by <?sin# (— gcosQ). The integration method we employ 
is a constant-timestep fifth order predictor-corrector [^9| . 

Fig. 1 shows a schematic drawing of the ball on the 
line. The spacing between two balls fixed on the line is 
2er, where e is a number which in the disordered case is 
chosen uniformly distributed in the interval [0, e max ]. The 
impact angle 7 is defined as the angle enclosed by the line 
joining the centers of the impacting balls and the normal 
to the plane. It is taken to be negative when the ball 
collides with the uphill facing side of a ball on the plane, 
positive on the downhill side. 
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FIG. 1. Schematic drawing of the ball on the line. 



III. SIMULATION RESULTS 

A. Global characteristics of the motion 

In all simulations presented here, the ball was launched 
onto the line with a rather high velocity v x in the x- 
direction and quite low v y . If inclination angle 9 and size 
ratio $ = R/r are in a suitable range, the moving ball 
very quickly, usually after passing only a few balls on the 
line, reaches a steady state with well-defined mean veloc- 
ities in the x- and y-directions. Clearly, the average over 
v v is zero in our problem, so the only interesting mean 
velocity is the average over v x , which we denote by v. To 
obtain this mean velocity, we first average over a certain 
number of timesteps (usually 500). This value is large 



enough to average out the comparably large fluctuations 
occurring during collisions, while it is still so small that 
the ball moves only a very short distance (much less than 
the radius of balls on the line) during this time. This av- 
eraged value clearly still gives a fluctuating velocity, but 
the fluctuations are very small, and v, the mean value of 
these averaged velocities, is well-defined. 

We investigated the motion of the ball both on lines 
with equally spaced balls and on lines with randomly 
spaced balls. We found that the essential features of the 
motion are alike in both cases. Fig. 2 shows the velocity 
of the rolling ball for various $ and equal spacing of balls 
on the line with e = as obtained from simulations using 
the above mentioned parameters. For the case e = 0, 
experimental data is available (though only for size ra- 
tios $ < 2.0) jl3|,[l6]]. Our stable mean velocities are of 
the same order of magnitude, but in the simulations the 
range of inclination angles for which v is well defined, 
i.e. a steady state is reached, is a bit narrower than in 
the experiments. This might be due to the difference 
in the experimental setup. There, a ball moved down a 
line of balls sitting in a V-shaped groove, and it is to be 
expected that contact with the groove walls influenced 
the motion as an additional source of dissipation due to 
friction and collisions with the walls, but it is unclear 
how strong this influence was. For a direct comparison 
of simulation and experiment, it would be desirable that 
a "more" 2-dimensional experiment be done, like for ex- 
ample a ball moving down a row of cylinders, to rule out 
these boundary effects. 




0.20 



FIG. 2. Dependence of the mean velocity v on the incli- 
nation of the line for various size ratios <&: $ = 1.5 (circles), 
$ = 2.0 (squares), $ = 2.5 (triangles), $ = 3.0 (stars). 

For angles lower than the one for which the smallest 
steady state velocity is reached for a given the ball 
loses all the initial velocity it had and very quickly stops, 
usually after passing only very few balls on the line (re- 
gion A of the phase diagram). We will denote this min- 
imum angle, for which a steady state with v ^ still 
exists, by 9 AB ($). This defines the phase boundary be- 
tween region A and B. 

All velocity curves in Fig. 2 exhibit a sudden increase 
of v to a value where it remains roughly constant. This 
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sudden increase has also been observed in experiments 
|p3| . For angles smaller than the one where this hap- 
pens, the velocity of the ball shows only very small fluc- 
tuations and the steady state velocity v does not depend 
on the initial velocity of the ball. At the inclination an- 
gle where the velocity suddenly shoots up, the behaviour 
changes qualitatively. The fluctuations of v x increase sig- 
nificantly, and the behaviour of the ball now depends on 
the initial velocity. Depending on the starting velocity, 
the ball can accelerate (usually if the starting velocity is 
larger than the stable mean velocity) and start to jump 
visibly, or it can reach a constant velocity (if released 
with an initial velocity smaller than the stable v in this 
region). However, even when the ball does not acceler- 
ate, in this region v can depend somewhat on the initial 
velocity. We thus denote the angle at which the sudden 
increase takes place by 9 BC (&), since it defines the upper 
boundary of the constant velocity region B independently 
of the initial velocity of the ball. 

If 9 increases even more, the ball accelerates and starts 
to jump significantly (the length and height of the jumps 
reaches a few ball diameters). We did not investigate 
this jumping motion any further, for the following two 
reasons. In the simulations the ball accelerated up to 50 
m/s and more, which is a velocity the ball would never 
reach in experiments due to air resistance. Thus, the 
question of whether the ball in the jumping regime can 
reach a steady state only due to collisions with the plane 
will not be discussed here, even though it is of theoretical 
interest. The second reason is that for grazing impacts 
at high velocities, an artefact of constant timestep algo- 
rithms, the so-called brake failure effect J2{| may set in, 
which leads to anomalous dissipation of energy. Its on- 
set can be shifted to higher velocities by increasing the 
spring constant k n , i.e. the stiffness of the balls. This, 
however, decreases the timestep, thus increasing the sim- 
ulation time tremendously. The constant timestep algo- 
rithm then becomes a very ineffective way of simulating 
the motion, since the ball spends most of its time in free 
flight, where a timestep small enough to integrate colli- 
sions correctly and to avoid spurious effects is essentially 
a waste of computation time. In this regime, event driven 
(ED) simulations would be more appropriate, which is 
why we will not discuss the motion in the "high bounce" 
regime. One might argue that, since the ball moves down 
the line in a series of bounces, an ED algorithm might 
in any case be a more efficient and appropriate way to 
simulate the motion. As we will see later, this is not the 
case, as both short and long-lasting contacts occur, the 
latter of which are not treatable by an ED algorithm in 
a straightforward way |2g] . 

The velocity curves obtained from experiments as well 
as from simulations suggest the functional form 



= «o($) + /($)Vsin0-sin0 AB ($) 



(8) 



on <f>. /($) denotes a (still unknown) scaling function, 
which, unlike in 3D, does not seem to be a simple power 
law. We thus plot (v — v Q (&)) 2 in Fig. 3, where i>o($) 
is obtained by fitting a square root to w(sin(9). vq(<&) 
typically is of the order of 4 cm/s. The error bars give 
the variance of the averaged velocities, averaged (as the 
plotted value of v itself) over a number of simulation runs 
with different starting velocities. Fig. 3a shows velocities 
in the case e = for various size ratios, Fig. 3b illustrates 
the influence of disorder on the velocity of a ball of size 
ratio $ = 2.25. 
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for a fixed value of r in region B. All curves start at a 
certain offset velocity vq, which seems to depend slightly 



FIG. 3. Dependence of the mean velocity v on the incli- 
nation of the line (a) for the same parameters as Fig. 2; (b) 
for various spacings of the balls on the line (<E> = 2.25): balls 
equally spaced with e = (stars), e = 0.1 (triangles), e = 0.2 
(circles); balls disordered with e max = 0.2 (squares). 

Fig. 3b demonstrates how the velocity of the ball is 
affected by the introduction of disorder to the line. It 
shows that v on the disordered line with e max = 0.2 can 
be approximated by the velocity on a line with equally 
spaced balls and a spacing e = 0.1, corresponding to the 
mean value of the disordered case. It also shows that 
BC ($) depends on both $ and the arrangement of the 
balls on the line, whereas the maximum v only seems to 
depend on $. 

In Fig. 4 we plot the corresponding phase diagram for 
three cases: two cases for a line with balls equally spaced, 
with e = and e = 0.2, the third for a disordered line 
with e max = 0.2. The lines denote the phase boundaries 
given by the angles 9 AB (&) and 9 BC (&) defined previously. 
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Obviously, the introduction of disorder has the same ef- 
fect on the phase boundary AB as the introduction of 
an equal spacing of balls with e = e max . This is under- 
standable; the stopping of a ball should be ruled by the 
deepest "traps" into which it might fall. The boundary 
BC rather seems to be determined by the mean spacing 
of balls, since in the disordered case it falls somewhere 
between the two extreme cases of an ordered array with 
e = and e = e max . 




FIG. 4. Phase diagram for e„ = 0.7 and various spac- 
ings of the balls on the line: balls equally spaced with e = 
(solid line) and e = 0.2 (dashed line); balls disordered with 
£max = 0.2 (dot-dashed line). 



B. Detailed dynamics of the motion 



and the corresponding impact velocities. The times of 
flight between impacts and corresponding previous im- 
pact angles are equally strongly correlated (see Fig. 6). 
From both Fig. 5 and 6 it is obvious that the ball moves 
down the line in a series of bounces. There is even a 
certain range of impact angles that are never hit. In ad- 
dition, we can extract from Fig. 6 that the typical times 
of flight between impacts correspond to a distance of less 
than 3 mm, which is smaller than the radius of the balls 
constituting the line. So the bounces the ball undergoes 
cannot be very high or far, and the bouncing ball collides 
with every ball on the line several times. 
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FIG. 5. Correlation of impact velocities and impact angles 
for $ = 4, sint9 = 0.05, e = 0. The upper points correspond 
to the relative translational tangential velocity v t = v ■ s, the 
lower ones to the relative normal velocity v n . 



In order to investigate the mechanism by which the 
ball maintains a constant velocity in region B, we have 
to look into the details of the motion. Certainly, distri- 
butions of the impact angle 7, of times of flight between 
impacts or of the impact velocity v 1 would be of interest. 
One might ask as well if there are correlations between 
impacts at certain angles and the corresponding impact 
velocities or between the time of flight after an impact 
and the corresponding impact angle. Actually, the latter 
provides even more detailed information; for this reason 
we will first have a look at these correlations. We will 
first discuss the ordered case of a line with no spacing 
between the balls, i.e. with e = 0. 

In Fig. 5 we plot the velocity of the ball right before 
an impact as a function of the corresponding impact an- 
gle 7. Each dot corresponds to a collision. We chose to 
plot the normal velocity v n and the tangential transla- 
tional velocity v t — v ■ s rather than v x and v y , as they 
provide more information on the mechanisms involved. 
Since dissipation takes place only through the normal 
velocity (dissipation due to F s is negligible, as we will 
see later) , the evolution of this quantity may explain the 
mechanism of energy loss maintaining the steady state. 
Since v t is the velocity component of the motion of the 
center of mass of the ball along the bumps of the surface, 
it reflects how this bumpyness is felt. It is obvious that 
there is a strong correlation between the impact angles 
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FIG. 6. Correlation of time of flight after an impact at an- 
gle 7 and corresponding previous impact angles (parameters 
as in Fig. 5). 

How can these results be understood? From Fig. 5 two 
important points can be extracted. Firstly we see that 
at certain impact angles 7, the moving ball is likely to 
hit the ball on the line with a well-defined corresponding 
normal velocity v n , which leads to an equally well-defined 
time of flight (see Fig. 6). Secondly, the normal velocities 
are largest for negative impact angles (i.e. hitting on the 
uphill side of a line ball) and get smaller (and eventually 
very close to zero) with increasing impact angles. The 
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second observation helps to explain the correlations and 
reveals the reasons for the regularity of the motion. 

Consider a ball that has just arrived at the maximum 
possible positive impact angle 7 max on a certain ball on 
the line, say ball number k. The value of 7 max is given by 
the geometry and defined as 



7- 



1 



arcsm ■ 



1 + $ 



(9) 



From Fig. 5 we see that at 7 max the moving ball has lost 
nearly all normal velocity with respect to ball k, i.e. it 
is rolling or sliding down the lower part of the "downhill 
side" of the corresponding ball on the line. Even though 
the total shear velocity v s at the point of contact is quite 
small (the ball rolls, thus rotational velocity is opposed to 
the translational tangential velocity and nearly compen- 
sates it), the translational tangential velocity v t — v ■ s, 
which we plot in Fig. 5, is quite large (close to v for 
larger $). Immediately after having reached 7 max on ball 
k, the moving ball impacts the next ball, k + 1, at — 7 max , 
where a large part of this previously tangential velocity 
is now normal velocity, as the direction of the vector n 
connecting the centers of the impacting balls with respect 
to v has changed. The ball thus gets thrown up again, 
and manages to reach the downhill side of ball k + 1 
after a few jumps. The comparably high normal veloc- 
ity of the ball at — 7 max is lost in two ways in crossing 
a fixed ball: most is lost by dissipation due to impacts, 
but partly it is converted into tangential velocity by the 
increasing obliqueness of successive impacts at positive 
7. In the process of bouncing over the top of ball k + 1 
on the line, the moving ball loses nearly all of its normal 
velocity, so that it again reaches 7 max with nearly only 
tangential velocity. These steps repeat themselves over 
and over again while the ball moves down the plane, thus 
retaining a constant mean velocity. Note that strong ge- 
ometrical constraints prevent the ball from rolling down 
the line without ever bouncing. Whenever the moving 
ball rolls down the downhill side of a ball k, it is thrown 
onto the uphill side of ball k + 1 with considerable normal 
velocity with respect to ball k+1 as a result of the change 
in the geometry. To remain in contact with ball k + 1, 
the moving ball would to have to lose a very substantial 
amount of this normal velocity in a single impact, or it 
would jump up. Persistent rolling thus is only possible 
in the case of vanishing normal restitution. Rolling on 
part of each ball is possible, however, as we will see in a 
moment. 

In order to describe things more quantitatively, we take 
a look at the distribution of impact angles. Fig. 7a gives 
an example for an angle of inclination in the middle of 
region B of the phase diagram for <E> = 2.5. The distri- 
bution exhibits clear peaks for 7 < 0. As 7 increases, 
the peaks broaden somewhat and approach each other 
until they are nearly indistinguishiblc in the histogram, 
although there still is structure in the correlation plots. 
For plane inclinations closer to 9 AB than 9 BCl the distri- 
bution, like in Fig. 7a, even breaks off at a value 7 < 7 max , 



indicating the start of a long-lasting contact. Here, the 
ball has lost so much of its normal velocity in impacts it 
suffered on crossing over the top of a ball on the line, that 
it finally starts to roll over part of this ball. This rolling 
motion can even start while the moving ball is still on the 
uphill facing side of the fixed ball, though this only takes 
place very close to 9 AB . All these impacts end at 7 max , 
which leads to the very pronounced peak at — 7 max . Inte- 
grating the distribution over negative and positive angles 
respectively yields the result that there are actually more 
impacts for positive than for negative 7. 
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FIG. 7. Distribution of impact angles for $ = 2.5, e = 
and (a) sin (9 = 0.057 (b) sin<9 = 0.075 (c) sin<9 = 0.091 (d) 
sine = 0.103. 

When 9 is increased, the peaks for negative 7, except 
the one at — 7 max , move towards zero, finally disappear- 
ing into the continuous distribution for positive 7 (sec 
Fig. 8). At a large enough inclination, they get visi- 
bly broader, showing that the velocity at 7 max is now no 
longer as sharply defined as before, probably due to the 
fact that not all the normal velocity can be lost on only 
one ball. Still, the ball does not accelerate as long as 
its starting velocity is only a little higher than its stable 
mean velocity (this is the region beyond the transition 
from linear to more irregular behaviour in Fig. 2). 

Fig. 7(b-d) shows how the motion of the ball changes 
qualitatively at 9 BC . The 9- values are taken a little be- 
low (b), as close as possible to (c) and a little above 9 BC 
(d). The heretofor clearly defined peak at — 7 max broad- 
ens considerably and even seems to develop a small side 
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peak as 6 reaches BC . The qualitative change of the be- 
haviour of the ball can also be observed in the velocities. 
In particular in the angle region where the velocity curve 
flattens out again, an intermittent behaviour of the ball 
can be observed - it will accelerate a little, even start to 
jump a little, but then suddenly get braked again. The 
mean velocity in this case is determined by how fast the 
ball is accelerated and braked respectively. The same in- 
termittent behaviour of the ball shortly before passing 
from motion with a mean constant velocity to a jumping 
regime has been observed in 3D experiments ]l6[ ] and in 
a stochastic model for the 2D case uM. 
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FIG. 8. Location of the peaks in the distribution of impact 
angles for $ = 2.25. The error bars denote the width of the 
peak. The solid, dotted, dashed, long-dashed and dot-dashed 
lines denote respectively the 1st, 2nd,..., 5th distinguishable 
peak, excluding the one at — 7 max . 



While the intermittent behaviour marks the transition 
from region B to region C in the phase diagram, stopping 
of the ball takes place when the typical tangential veloc- 
ity at 7 max does not convert into enough normal velocity 
at — 7 max to carry the ball over the top of the next ball in 
the line in a few jumps, or at least up to a point where 
the remaining tangential velocity suffices to make it roll 
over the top of this ball. 

In regime B, the motion is not only characterized by 
typical impact angles 7, but also by a very strong cor- 
relation of successive impact angles. In Fig. 9, we plot 
7„+i, the impact angle for impact n + 1 as a function of 
the previous impact angle 7„. Even when introducing a 
random spacing of balls on the line the strong correla- 
tion remains. This observation leads us to the discussion 
of the behaviour of the moving ball on a disordered line. 
The phase diagram already suggests that introducing dis- 
order has a similar influence on the motion of the ball as 
introduction of a regular spacing. In both cases, the re- 
gion of stable motion in the phase diagram shifts to larger 
$ and 9. Plotting the distribution of impact angles for 
the cases displayed in Fig. 3b shows the effect of disorder 
(see Fig. 10). 



FIG. 9. Correlation of successive impact angles for $ = 
2.25, sinf? = 0.094. (a) Balls on the line regularly spaced 
(e = 0), (b) balls on the line disordered with e max = 0.2. 



We checked that this regularity in the case of a dis- 
ordered line is not a finite size effect. The angle dis- 
tributions do not change when the size of the system is 
increased. From Fig. 10 it can be seen that the distri- 
butions for the impact angles in the disordered case lie 
between the limiting cases e = and e = 0.2 of an or- 
dered line. In addition, the peaks for e = 0.1 correspond 
to the center of the peaks in the disordered case. The 
reason for the large width of the peaks in the case e = 
lies in the fact that here the velocity is already quite 
high (see Fig. 3b). We interpret the results of the disor- 
dered case in the following way. The motion of the ball 
is influenced mainly by two factors - the minimum and 
maximum spacing of balls on the line. If the ball is to 
keep its mean velocity, without being stopped and with- 
out being accelerated, the velocity has to have a value 
that enables it to get out of the "deepest valleys" existing 
between two balls (given by e max ), but still is low enough 
that in all cases (even in those where e = 0) most of the 
normal velocity is dissipated on crossing over the corre- 
sponding ball so that at 7 max of this ball only tangential 
velocity is left. The variations in this tangential velocity 
are so small that they only broaden the peaks for the 
impact angles, but do not lead to a qualitative change of 
the motion of the ball. 
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FIG. 10. Distribution of impact angles for $ = 2.25, 
sin# = 0.088 and various spacings of the balls on the line, 
(a) regular spacing, e — (b) disordered line, e max = 0.2 (c) 
regular spacing, e = 0.2 (d) regular spacing, e — 0.1. 



C. Influence of material properties on the motion 

Having uncovered the mechanism by which the ball 
moves down the line and keeps a constant average ve- 
locity, we now have to ask how much this mechanism 
and the global results like v are influenced by material 
properties. The material properties incorporated in our 
simulations are the coefficient of normal restitution e n 
and the friction coefficient /i. We will show in this sub- 
section that their influence on v and the mean rotational 
velocity Co is very small. 

Experiments in 3D already indicate that the charac- 
teristics of the motion are hardly influenced by material 
properties |L4]]27]]. Rolling steel, glass and plastic balls 
down the plane gives nearly the same velocity, though 
for plastic (which has the lowest e„ and largest fj) the 
B region of the phase diagram is found to be somewhat 
extended. Our simulations in 2D indeed show that the 
mean velocity v is nearly independent of both e„ and fi. 
Fig. 11 demonstrates this for the case of a ball of size 
ratio <& = 2.25. Fig. 11a shows the velocities for vary- 
ing normal coefficient of restitution e n . It can clearly be 
seen that e„ influences the phase diagram, i.e. the exten- 
sion of region B, but has only a very small influence on 
v. Though 9 AB is hardly affected by e„ (except for very 
small $), 6> BC moves to larger inclination angles with 
decreasing e„. But e„ influences v slightly in a direction 



that is contrary to what one would expect intuitively. In- 
creasing the dissipation leads to a slight increase in the 
velocity. We have found no explanation for this so far. 
We have, however, understood the relative insensitivity 
of v to e n . 
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FIG. 11. (a) Dependence of v on the dissipation for 
$ = 2.25 and e = 0: e„ = 0.4(O), 0.5(a), 0.6(O), 0.8(A). 
(b) Dependence of v and luR on the friction coefficient y,: fi = 
0.1(O), 0.3(D), 0.5(O), 0.7(A), 1.0(*). 

The reason for the small influence of e n is essentially 
the regularity of the motion. As shown in the previous 
subsection, the ball moves over each line ball in a suc- 
cession of bounces, in the course of which it loses all or 
most of its relative normal velocity with respect to this 
line ball. e„ mainly determines how many bounces are 
necessary to achieve this. As long as the moving ball 
has only a negligible amount of normal velocity left as it 
reaches 7 max , it seems to be unimportant for the tangen- 
tial velocity at this point how many impacts were needed 
to lose the normal velocity the ball had at — 7 max . 

For the phase boundary BC, the value of e n is impor- 
tant. Since we found that the motion of the ball starts 
to get unstable when the motion is no longer regular, we 
expect the destabilization to start when the normal ve- 
locity at — 7„ lax with respect to one ball cannot be lost by 
the time 7 max is reached. But e„, which determines the 
energy loss in a collision and thus the height (and thereby 
the length) of the next jump, together with v determines 
how many jumps are made on one line ball. It thus gives 
an upper limit to the maximum amount of energy that 
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can be dissipated on a single line ball. The smaller e„, 
the more energy is dissipated in each jump, and the more 
jumps are possible, as their length decreases. Thus, 9 BC 
shifts to larger 9 for a given $ with decreasing e„ and 
approaches the static angle of stability, which in the 2D 
case is given by 7 max . 

We find that because the ball is allowed to rotate, the 
value of the friction coefficient /i is even less important for 
the behaviour of the ball than e„ (see Fig. lib), /i influ- 
ences neither v nor the mean rotational velocity u> signifi- 
cantly. It also does not change the phase boundaries. The 
reason for this is that for the range of 9 considered here, 
the ball can be expected to roll without slipping most of 
the time independent of ^ if our implementation of the 
tangential force law correctly reproduces Coulomb's law 
for sliding friction. Let us assume for a moment that 
the ball, while moving down the line, is always in con- 
tact with the balls on this line, which is only possible if 
e n = 0. It would then roll without slipping if at all 7 the 
condition \F S \ < n\F n \ were fulfilled. For a rolling sphere 
under the action of gravity pj| , 

\F s \ = l\g-A, (10) 



IV. A THEORETICAL MODEL FOR THE MEAN 
VELOCITY 

In this section, we derive v by a simple analytical treat- 
ment. Our simulation results play an important role 
here, as they show which simplifications can be intro- 
duced without losing essential features of the motion. 

From our results on the influence of the coefficient of 
restitution e„ on the motion of the ball, we can deduce 
a very accurate result for v in the case of equally spaced 
balls on the line. Since we find that v hardly depends 
on e n we investigate the limiting case e n = 0. As we 
also found that, due to friction, the ball is able to sus- 
tain v s = on average (i.e. the ball would roll without 
slipping in the case e n = for large fi), we make the 
following assumptions. 

We assume that the moving ball is always in contact 
with the fixed balls, i.e. rolls down the line without slip- 
ping. Thus at all times, 

|«1 =v t =v-s, (12) 

V n = 0, (13) 



so that the criterion for fi for the ball to roll without 
slipping reads 

M>^|tan(0 + 7 )|. (11) 

Though close to 7 max this condition is usually not fulfilled 
(the smaller 4>, the larger the region of 7 for which slip 
can occur), it holds on the largest part of a ball on the 
line even for small fi. We would thus assume the ball on 
the average to roll without slipping in the case of e n = 
even for small values of the friction coefficient fi. 

In our simulations, however, we used larger values for 
e n , like e n = 0.7 in Fig. lib, so there the ball rather 
bounced than rolled down the line. But the distances 
the ball covers between bounces in the steady state are 
very small. We find in our simulations that the ball, 
which we launch onto the plane without rotational ve- 
locity, soon picks up rotation during impacts, such that 
when the steady state is reached, the rotational velocity 
has adjusted itself to a value which on average leads to 
zero relative velocity of the surfaces of the moving ball 
and the fixed balls (except close to 7 max , just as would 
be expected from eq. (pi])). In the free flight between 
collisions, the ball picks up translational velocity while 
the rotational velocity remains unchanged. The moving 
ball has thus gained excess shear velocity, which is con- 
verted nearly completely into rotation in the next impact, 
if it was small enough. This is usually the case in the 
steady state, where distances covered between bounces 
are small. 



uiR = v t . (14) 

From these conditions the equation of motion of the 
ball on a single ball on the line from — 7 max to 7 max can 
be derived. The kinetic energy of the ball is 

£ kln = X -mv 2 t + X -Ju? (15) 

where J denotes the moment of inertia of the moving 
particle, which in the case of a sphere of radius R takes 
the value J = |mi? 2 . Using condition (^J), we get 

E kin = ^mv 2 t (16) 

for a sphere. The potential energy depends on the loca- 
tion of the moving ball on the line ball: 

E pot =m 5 (i? + r)cos( 7 + #). (17) 

Since vt = (R + r)j, the energy balance reads (denoting 
the energies at the start of the motion by E®. n and E® ot ) 

EL - ^rn(R + r) 2 7 2 = mg(R + r) cos( 7 + 9) - E° pot . 

(18) 

Differentiating with respect to time yields the equation 
of motion for j(t) 

7 = y-^sm( 7 + 0). (19) 

The completely inelastic collision at — 7 max that occurs 
when the moving ball passes from one line ball to the next 
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defines the boundary conditions for the problem. We de- 
note the velocity of the ball at 7 max by «t(7 max ) = w;. 
Since we assume rolling without slipping, on this ball on 
the line we had 



LOiR 



(20) 



with uji = w(7 max ). In the next instant, the moving ball 
hits the next ball on the line at — 7 max , with the tangen- 
tial velocity vt{— 7 max ), which with respect to this ball is 
(due to the change in geometry) 



7mw) = Vi C0S(27max). 



(21) 



In the collision that is now about to take place, the nor- 
mal component 



>(-7max) = WiSin(27max) 



(22) 



is reduced to zero, but this is not the whole effect of the 
collision. Since v t dropped before the impact due to the 
change in geometry, but ui was left unaffected, at — 7„ lax , 
there is excess rotational velocity and thus excess shear 
velocity v s . If the frictional force is high enough (which 
we assume in this treatment), then the shear velocity 
at the point of contact should again be reduced to zero 
during the impact, such that tangential and rotational 
velocity after the impact at — 7 max , which we will denote 
by vt and u>f respectively, fulfill the condition 



(23) 



During the impact at — 7 max the rotational and transla- 
tional velocities thus adjust themselves, with only negli- 
gible energy loss (friction here mainly helps to distribute 
the excess rotational velocity to the translational degree 
of freedom, but, as the ball can rotate freely, dissipates 
only very little energy during the adjustment). Before 
the adjustment of rotational and translational velocities, 
the kinetic energy of the system was 
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E Un = -mvf cos 2 (2 7max ) + -Jw. 



(24) 



With eqs. (|16|) and (J23J) we thus obtain from the energy 
balance the condition for Vf 
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-i; 2 cos 2 (2 7max ) + -Jujf = y-mt$. 

Substituting w$ according to eq. (Boh, we get 



v f =^/-(5cos 2 (27 max ) + 2) 



(25) 



(26) 



for the tangential velocity after the collision at — 7 max . 
Since v t = (R + r)j, this provides the boundary condi- 
tion for eq. (pj|). 

We obtain v numerically by starting at — 7 max with an 
arbitrary and very high starting velocity 70, letting the 
system evolve according to eq. (|T^). Whenever — 7 max 



is reached, eq. (|26J) is applied, and we start again at 
7max- We continue this process until cither the ball rolls 
back (which we consider as trapping) or until it reaches 
a steady state. The mean velocity in this steady state is 
plotted in Fig. 12, in comparison with simulation data for 
various coefficients of restitution. Clearly, the simulation 
for small e„ is closest to the theoretical curve (which as- 
sumes e„ = 0), but also for higher e n the approximation 
is still good. 

For the case e > eq. ( |l9| ) and (26) yield equally good 
results as in the case e = presented in Fig. 12. In the 
case of a disordered line, an estimation of v with the value 
of 7 max chosen according to the mean value of e fits the 
simulation results equally well. 




FIG. 12. Theoretical prediction of v (solid lines) for e = 
under assumption of e n = and simulation data for e n = 0.1 
and size ratios $ =1.75 (circles), 2.25 (squares), 3.0 (trian- 
gles) . The dashed and dot-dashed lines respectively also show 
the corresponding simulation data for e n = 0.5 and 0.7. 

Another result that can be obtained from this theoret- 
ical treatment is the phase boundary 9 AB . To this end, 
we make use of eq. (Q). By setting 



and 



E 



7 

To 



m(i? + r) 2 7 2 (- 7max ) 



mg(R 



7»ax), 



(27) 



(28) 



the values at the begin of the motion over one ball in the 
line, and by using the fact that in the steady state 



7(-7max) = e t 7(7 max ) 



with 



(5cos 2 (2 7max ) + 2), 



(29) 



(30) 



we obtain from eq. 
the steady state 

^(^ax) 



the starting velocity j(— 7 max ) in 
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7 R + r 1 



~2 sin 7 max sin ( 



(31) 
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We now assume that the phase boundary 9 AB is reached 
when the moving ball arrives at the angle 7 = — 9 with 
zero velocity, since from there it can roll down simply 
by the action of gravity, even with zero starting veloc- 
ity at this point. For smaller inclination angles, the ball 
would roll back before reaching this point and thus stop, 
for larger inclination angles it would pass this point with 
some velocity and move on. From eq. ([l^) we obtain 
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7 R 



(1 — cos( 



7—)) (32) 



by setting 7 = — 9 AB and j(— 9 AB ) = 0. If the ball is in 
region B, the steady state condition ( |3l| ) has to be ful- 
filled as well, so that from eq. ( |3l| ) and eq. ( |32| ) we get 
an equation for 9 AB 



2e 2 



■ sin 7„ 



sin 9 A 



which finally yields 



sin A 



sin 7 maj 



. 2 

sin 7« 



(33) 



(34) 



cos 2 7„ 



1-e 2 



2e t 



1 



— „ cos 7„ 



Fig. 13 shows a comparison of this result to simulation 
data. Obviously, our theoretical result approximates the 
simulation results best for larger $. For <3? close to 1 the 
deviations get quite large, since here the assumption that 
the ball loses its normal velocity in a single impact is not 
fulfilled any more as well as before even for e„ = 0.1. In 
any case, eq. (poj) provides a lower limit for the value of 
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FIG. 13. The phase boundary 8 AB for e n = 0.7 (dotted 
line), e n =0.1 (dashed line) and the theoretical result from 
eq. ( piil ) for e n = (solid line). 

After the completion of this work, we became aware 
of work by Ancey et al. |2{| paralleling our theoretical 
treatment. They, however, neglect the effect of the rota- 
tional velocity in the impact at — 7 max , though rotation 
is explicitly included in their equation of motion. They 



thus have to introduce a fitting parameter to match their 
curves to experimental data. This is not necessary if the 
influence of rotation is included in the boundary condi- 
tions of the problem, as we have shown. We emphasize 
that eqs. ( |l9| ) and ( prjj ) hold for the instantaneous veloc- 
ity only in the case of vanishing normal restitution and 
perfectly rough surfaces, conditions which are not likely 
to be fulfilled by any commonly used material. However, 
as we have shown, the mean velocity of the ball is not 
influenced by the coefficients of restitution and friction 
and thus is correctly described. 

One more result can be deduced directly from our 
theoretical treatment. Since 7 ~ (r($ + I))- 1 / 2 , and 
v t = ($ + 1)^7, we expect v scale with t/t, i.e. to de- 
pend on the absolute size of the balls. Note that it does 
not, however, simply scale with \/W+l as well, as $ also 
enters into the boundary conditions and thus changes 7. 
But, keeping the geometry of the system (i.e. <!> and e) 
constant, we expect a scaling of v with i/r in our simu- 
lations. 

This scaling with r is observed in our simulations (see 
Fig. 14). In 2D experiments, it had been observed in 
an experimental setup restricted to $ = 1 fl5|j , here we 
find it as well if we keep $ fixed at an arbitrary value. 
This scaling seems to be quite universal, since it was also 
found in 3D experiments |30| and in the 2D stochastic 
model fill. 
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FIG. 14. Scaling of v with r keeping $ constant, r =0.5 
mm (o), 1.25 mm (A), 5 mm (Q). 



V. SUMMARY AND DISCUSSION 

We have discussed the mechanism stabilizing the mo- 
tion of a ball on an inclined line consisting of equally sized 
balls. Two important results emerged from our simula- 
tions. Firstly, we found that the motion of the ball in 
the steady state is very regular and consists of a series 
of small bounces on each ball on the line, contrary to 
what so far has been assumed Jl^Jl^]. In the course of 
these bounces the moving ball loses all relative normal 
velocity with respect to this ball. This dissipation mech- 
anism holds both for lines with equally spaced balls and 
for lines with a random spacing of balls. Furthermore, 
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the location of these bounces is the same on each ball, so 
that a random spacing of the balls on the line only smears 
out these locations a bit, but does not alter the motion 
significantly. Thus, the velocity of the ball in the case 
of a random spacing of balls on the line can be approxi- 
mated by the motion of the ball on an ordered line with 
appropriate spacing. The main reason for this regularity 
of the motion is that the moving ball has to climb over 
the top of every ball forming the line by a few bounces 
in the steady state. Since disorder only slightly changes 
the height to be climbed, it has only a small influence on 
the motion. 

Clearly, this mechanism for keeping a constant velocity 
cannot hold in 3D, where there is a number of ways for 
the ball to choose to go down the plane. The 3D case 
should be ruled by a competition between stability prop- 
erties of the plane and the preferential direction given by 
the plane inclination and inertia of the moving ball. One 
might argue that introducing stronger disorder in the 2D 
case, for example by using polydisperse balls on the line, 
would model a longitudinal section of the plane much 
closer. But this will not eliminate the strong dimensional 
difference between the 2D and 3D cases. In 3D, the mov- 
ing ball, hitting a ball on the plane a little on the side 
(which would be one of the smaller balls in the 2D sec- 
tion), would be deflected towards the side, which in 2D 
is impossible. In the constant velocity regime in 2D on 
a line with size polydispersity, the moving ball will again 
lock into some kind of quasiperiodic behaviour for the 
same reasons as explained above. The problem of hav- 
ing to overcome some maximum threshold (to get out of 
the deepest valley between two balls on the line) to keep 
moving is even stronger, contrary to the 3D case. There a 
ball, having suffered a large impact which brakes it a lot, 
might still find some way around the bump that braked 
it and accelerate until the next impact via some stati- 
cally unstable path. Thus, the motion in 3D is expected 
to be far more irregular than can be modeled in 2D. Fu- 
ture simulations in 3D will have to show where these 
differences originate and what the stabilizing mechanism 
consists of in this case. A stochastic model for the 2D 
case Jl8|,^4| nevertheless showed a viscous friction force, 
i.e. a linear dependence of v on sin0. This might be due 
to the fact that the randomness introduced there in the 
choice of the next impact captures the randomness of the 
"real" 2-dimensional plane. However, this question can 
only be answered ultimately by 3D simulations. 

The second important result is that material proper- 
ties like the normal coefficient of restitution e„ and the 
Coulomb friction coefficient /j, hardly influence the mean 
velocity v in the steady state. From this result and the 
knowledge of the mechanism stabilizing the motion of the 
ball in 2D we were able to predict v theoretically under 
the assumption of completely inelastic collisions. Only 
geometrical considerations sufficed for this treatment and 
gave a very good approximation to our simulation results. 
The results in this extreme case enabled us to derive a 
lower limit for the phase boundary 6 AB , which for large 



<!> gives a good approximation of # AB regardless of the 
value of the coefficient of restitution. This result should 
be relevant to the problem of segregation in the flow on 
inclined planes or in rotating drums. 

Though we have made significant progress in under- 
standing the stabilizing mechanism by which the ball 
keeps its constant velocity in 2D, a few open questions 
remain. Due to the strong nonlinearity of the prob- 
lem (computation of successive impact angles, velocities 
etc. from the equations of motion would require the so- 
lution of a fourth order polynomial), a direct iteration of 
the equations of motions as for e n = is very cumber- 
some, though in principle possible, in the case of non- 
vanishing e n . Since the phase boundary 9 UC , which sep- 
arates steady motion of the ball from the accelerating 
regime, depends on e„, it cannot be derived directly from 
our theoretical treatment. The most obvious simplifica- 
tion one might introduce, namely ignoring the structure 
of the plane in the computation of the next impact, is 
out of the question, since we have found this structure 
to be essential for the stabilizing mechanism. The prob- 
lem seems to be related to the classical problem of a ball 
bouncing on a vibrating plate fjlf , which lately has been 
discussed for the case of a partially inelastic ball |3^] . In 
the case of finite restitution, it was found that neglecting 
the motion of the vibrating plate can lead to erraneous re- 
sults, like the observation of "chaos" in a region of phase 
space where a more exact treatment shows the existence 
of eventually periodic orbits. Unfortunately, even this 
simple one-dimensional problem can only be solved by 
approximations in the extreme cases e„«l and e„ — > 1 , 
so it is by no means obvious how to relate these results 
to our problem, since for these extreme cases, we already 
have solved the problem. (The case e„ — > 1 is trivial, 
since there region B of the phase diagram vanishes.) 

So for the 2D case, the remaining open questions are 
the following. So far, we do not understand well how 
the ball gets braked down to its stable velocity when it 
is launched on the plane with a much higher initial ve- 
locity, or why it accelerates as soon as it bounces only 
approximately once per ball on the line. In addition, 
the questions of whether in the bouncing region a ste ady 
state is reached as well, as has been suggested before Q , 
and what the characteristics of this steady state are, are 
of interest. Does the ball really move in a "chaotic" way 
in the bouncing regime, as is suspected |L3j, or might an 
eventually periodic motion exist as in the case of the ball 
on the vibrating plate, which due to the long transient 
so far could not be observed? 

Although many of the interesting features of granu- 
lar flow come about by the collective behaviour and the 
interaction of many particles, we have shown that even 
the elementary processes involving only a single particle 
moving in a dissipative, but fixed environment, can of- 
fer quite a few new insights and open up new questions. 
One of these questions, which is also of relevance to gran- 
ular flows in general, is that of dimensionality. It has al- 
ways been implicitly assumed that 2-dimensional model 
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systems are an adequate tool to study granular flows in 
general, i.e. that the qualitative behaviour will not be 
different in 3D. Our results show that even though this 
is probably true for behaviour in the bulk, where grain 
motion is very confined, care has to be taken in relating 
the behaviour on free surfaces in 2D to that in 3D. 
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